%% it computes the invariand distribution of b,Z
h_00     = zeros(Ntot,1);
inde_hig = max(2,find(b_vec_mat>b_up,1));
inde_low = max(2,find(b_vec_mat>b_dn,1));
inde_mid = max(1,find(b_vec_mat>b_md,1));

inde_z_up   = find(z_vec>=log(zhig),1);
inde_z_md   = find(z_vec>=log(zmid),1);
inde_z_dn   = find(z_vec>=log(zlow),1);

h_00((inde_hig-1)*N_z+inde_z_dn)  =q_hig*incid;
h_00((inde_low-1)*N_z+inde_z_up)  =q_low*incid;
h_00((inde_mid-1)*N_z+inde_z_md)  =q_mid*incid;

if incid<1
    inde_hig = find(b_vec_mat>b_up2,1);
    inde_low = find(b_vec_mat>b_dn2,1);
    inde_mid = find(b_vec_mat>b_md2,1);
    h_00((inde_hig-1)*N_z+inde_z_dn)  =h_00((inde_hig-1)*N_z+inde_z_dn)+q_hig*(1-incid);
    h_00((inde_low-1)*N_z+inde_z_up)  =h_00((inde_low-1)*N_z+inde_z_up)+q_low*(1-incid);
    h_00((inde_mid-1)*N_z+inde_z_md)  =h_00((inde_mid-1)*N_z+inde_z_md)+q_mid*(1-incid);
end

%% policy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if policy_expe==1 %subsidy experiment
    drift_I   =  max(0,fnval(cub_spl_d,[b_vec_mat;A*ones(1,length(b_vec_mat))]));
else % use steady state policy
    drift_I   =  max(0,ppval(pp_drift,b_vec_mat));
end

if inde_dilution ==1
    if policy_expe==1 %subsidy experiment
        ell        =fnval(cub_spl_l,[b_vec_mat;A*ones(1,length(b_vec_mat))]);
    else % use steady state policy
        ell        =ppval(pp_ell,b_vec_mat);
    end
else
    ell        = zeros(1,length(b_vec_mat));
end

inde_out          = length(find(b_vec_mat>b_bar));
inde_bb           = find(bb_vec_mat>b_bar);
inde_bi           = find(bb_vec_mat<=b_bar);
inde_nego         = find(b_vec_mat>=b_bar*a,1);

%% transition matrices  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inde_b            = reshape(ones(N_z,1)*(ell./b_vec_mat-(rho+drift_I-sigma^2/2)),Ntot,1);
Db                = Db_f;
Db(inde_b<=0,:)   = Db_b(inde_b<=0,:);
Db_db             = sparse(Db.*inde_b);
Dbb_dbb           = sparse(Dbb*(sigma^2/2));

inde_z            = reshape(ones(N_z,1)*((drift_I-sigma^2/2)),Ntot,1);
Tz                = Tz_f;
Tz(inde_z<=0,:)   = Tz_b(inde_z<=0,:);
inde_z_min        = 1:N_z:Ntot;
inde_z_max        = N_z:N_z:Ntot;
Tz(inde_z_min,:)  = 0;%impose that at minimum of grid f(zmin)=f(zmin+dz)
Tz(inde_z_max,:)  = 0;%impose that at maximum of grid f(zmax)=f(zmax+dz)
Tz_dz             = sparse(Tz.*inde_z);
Tzz_dzz           = sparse(Tzz*(sigma^2/2));

Tzb_dzb           = sparse(Dbz)*(sigma^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AA            = Db_db  + Tz_dz    + Dbb_dbb      - Tzb_dzb       + Tzz_dzz;
A_exit        = speye(Ntot)*delta;

%construct matrix to relocate upon renegotiation
E = zeros(N_z,Ntot);
for ii=1:N_z
    E(ii,(inde_nego-1)*N_z + ii)=1;
end
EE  = (repmat(sparse(E),N_b,1));%matrix of zeros, except for the location of the renegotiation

%construct the mass of exiters from each greed point j
X   = sum(AA(:,inde_bb),2);X(inde_bb)=0; 
XX  = repmat(X,1,Ntot);%repeat the vector in multiple equal columns

%obtain the transition matrix with renegotiation
AA2 = -A_exit + AA + (phi*(EE.*XX));%here add mass only when XX(j,i) is positive, which is when there is exit from point (j) AND EE(j,i)=1 that is when i is a renegotiation point given j
AA2(inde_bb,inde_bb)=-speye(length(inde_bb));

%transpose the transition matrix so that you go from HJB to KFE; the minus sign is to take it to left of the equation
A_T = -AA2';

%% find the distribution
inde_b_bar=find(b_vec_mat<=b_bar,1,'last');

crit_h=1;
if calib==1,sum_h=1;end
h_til     = zeros(Ntot,1);
while crit_h>1/100000
    h_til(inde_bi)  = A_T(inde_bi,inde_bi)\(h_00(inde_bi)*(1+sum_h*delta));
    crit_h =max(abs(sum(h_til)-sum_h));
    sum_h=sum(h_til);
end
f_dist = h_til;